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We study the low-energy dynamics of 5 = 1/2 antiferromagnetic Heisenberg clusters constructed 
by diluting a square lattice at vacancy concentration p at and below the percolation threshold 
p* ~ 0.407. The finite-size scaling behavior of the average excitation gap, (A) ~ , where 
L is the cluster length, is obtained using quantum Monte Carlo results for an upper bound A* 
to A, derived from sum rules. At the percolation threshold, we obtain a dynamic exponent z = 
3.6 ± 0.1 f» 2Df for clusters with singlet [S — 0) ground state. Here Df = 91/48 is the fractal 
dimensionality of the percolating cluster. We argue that this large dynamic exponent — roughly 
twice that expected for quantum-rotor excitations — is a consequence of weakly interacting localized 
effective magnetic moments, which form due to local sublattice imbalance. This picture is supported 
by an extremal-value analysis of local spectral gaps, which delivers an exponent relation (between 
z and two exponents characterizing the local gap distribution) reproduced by our simulation data. 
However, the average (A*) over all clusters, which have mostly ground state spin S > 0, scales with 
a smaller exponent than for the S — clusters alone; z ~ 1.5-D/. Lanczos exact diagonalization for 
small clusters show that typically, S S — I in the lowest-energy excitations, while the dominant 
spectral weight originates from S ^ S + 1 excitations. Thus, the scaling of (A*) for clusters 
with ground state S > does not reflect the lowest-energy excitations, but the higher S ^ S + 1 
excitations. This result can be understood within a valence-bond picture. To further explore the 
scenario of localized moments, we introduce a classical dimer-monomer aggregation model to study 
the distribution of nearest-neighbor sites forming dimers (which are the objects used in mapping 
to the quantum-rotor model) and unpaired spins (monomers). The monomers are localized, and, 
thus, effective magnetic moments should form in the spin system. We also study the lowest triplet 
excitation of S = clusters using quantum Monte Carlo calculations in the valence bond basis. 
The triplet is concentrated at some of the classical monomer regions, conflrming the mechanism of 
moment formation. The number of spins (and moment regions) affected by the excitation scales 
as a non-trivial power of the cluster size. For a dimer-diluted bilayer Heisenberg model with weak 
inter- layer coupling (where the system remains Neel ordered), there is no sublattice imbalance. In 
this case we find z Df, consistent with quantum rotor excitations. For a single layer at p < p* we 
find z « 2 = D, which indicates that the weakly interacting localized moment mechanism is valid 
only exactly at the percolation point. There is a cross-over behavior close to the percolation point. 

PACS numbers: 75.40.Gb, 75.10.Jm, 75.10.Nr, 75.40.Mg 



I. INTRODUCTION 



Two-dimensional (2D) antiferromagnets under dop- 
ing with non-magnetic impurities have emerged as in- 
teresting systems with rich possibilities to explore vari- 
ous disorder-driven phase transitions belonging to differ- 
ent universality classes J^'^'^'^'^ Non-magnetic impurities 
(vacancies) enhance quantum fluctuation by reducing the 
connectivity of the spins. Many earlier calculations^ for 
the 2D S — 1/2 Heisenberg model had indicated that 
the quantum fluctuation can become strong enough to 
destroy the antiferromagnetic long-range order at a va- 
cancy concentration pc less than the classical percola- 
tion threshold p* — whence Pc would be a quantum crit- 
ical point. However, more recent quantum Monte Carlo 
(QMC) simulations of the diluted quantum Heisenberg 
model^''^ studies of effective classical systemSfi as well 
as experiments on La2Cui_a;Znj;04 (with non-magnetic 
Zn substituting S" = 1/2 Cu ions)^" all suggest that long 
range order actually survives all the way up to the per- 
colation point p*, i.e., Pc = p* for the single 2D layer. 



The percolating cluster at p* is ordered^S which implies 
that the static properties at the dilution-driven transition 
in the quantum Heisenberg model scale as in the classical 
(percolation) problem. However, quantum fluctuations 
lead to changes in the low-energy spin dynamics. The 
critical exponents therefore in general depend on classical 
percolation exponents as well as the dynamic exponent 
z of the quantum spin clustersi^ The dynamic exponent 
of the percolating cluster is therefore important, and the 
focus of this paper. 

The dynamic exponent governs the scaling of the gap 
A between the ground state and the lowest excited state 
of a finite cluster. With L denoting the cluster length (de- 
fined in some suitable way for a random cluster with ir- 
regular shape), the gap scales, on average, as (A) ^ L^^. 
For a clean _D-dimensional antiferromagnetic system on 
a bipartite lattice with N (even) sites, every spin can 
be paired up with a nearest-neighbor spin on the op- 
posite sublattice to effectively form a "quantum rotor" 
with angular momentum / = 0, 1 states. In the map- 
ping to a quantum rotor model, these local degrees of 
freedom are replaced with angular momenta li taking all 
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integer values, with the high k states suppressed due to 
their energy being (x If. The ground state of the cou- 
pled quantum rotor system is a singlet. If the system is 
long-range ordered (but the global rotational symmetry 
has not been broken by any external perturbation) , then 
the low-energy excitations of the coupled rotors (and the 
Neel ordered spin system^^ ) are those of a single quantum 
rotor with mass oc N. Thus A ~ N~^, i.e., z = D. 

According to one recently proposed scenario for ran- 
domly diluted antiferromagnets,^ the quantum rotor 
states remain the lowest-energy excitations even at p*, 
where the dimensionality -D/ of the percolating cluster 
is fractal; z = Df = 91/48 Following the discussion 
above, this would seem to require that each spin can be 
paired up into a dimer with one of its nearest neighbors 
to effectively form a quantum rotor with I = ground 
state. This situation can be realized in the special case 
of the dimer-diluted bilayer,-'' in which two coupled layers 
are diluted exactly in the same way by removing inter- 
layer spin dimers. All the remaining spins can then be 
paired with spins on the opposite layer. At sufficiently 
weak inter-layer coupling, the ground state of the largest 
connected cluster of spins in this system is long-range or- 
dered for p < p* ,^ and, thus, the ground state should fall 
into the class of quantum rotor states with gap cx N^^. 
However, in the case of a single diluted layer (or a bi- 
layer with inter-layer coupling J± =0), there are in gen- 
eral some "dangling spins" (or more generally, regions 
with local sublattice imbalance) in which not all spins 
can be simultaneously paired up into nearest-neighbor 
dimers. One may still be able to pair spins over longer 
distances (which would also imply longer-range interac- 
tions between the rotors in the effective model), but at 
some point, when very long distances are required, the 
mapping to simple quantum rotors should break down. 

Our assertion is that, at the percolation point, there 
are regions of spins that effectively form isolated mag- 
netic moments, which cannot be described within an ef- 
fective model containing only coupled rotors. The spatial 
distribution of these moment regions, and weak effective 
interactions between them (mediated by the magnetically 
inert parts of the percolating cluster), lead low-energy ex- 
citations which are dramatically different from those of 
the quantum rotor system. We introduced this scenario 
and presented supporting numerical evidence in a recent 
paperji^ Using finite-size scaling, we found a consider- 
ably larger dynamic exponent than the quantum-rotor 
value; z « 2iD^ instead of z — Df. Here we provide more 
details of this work, and also expand significantly on the 
previous calculations. We use several different methods 
to indirectly and directly examine the low-energy exci- 
tations of different types of clusters, both at and away 
from the percolation point. 

The conclusion that z « 2Df for clusters at the per- 
colation point is based largely on quantum Monte Carlo 
(QMC) calculations of an upper-bound A* to the lowest 
excitation gap A for finite clusters with singlet {S = 0) 
ground states. The bound is defined using standard sum 



rules, discussed in detail in Sec. lHBl [and summarized as 
Eqs. ®, dZl), and ^]. The bound is exact, A* = A, 
for a spectrum with a single mode, and is known to scale 
with the system size in the same way as A more gen- 
erally, e.g., in the clean Heisenberg model^S. It can be 
evaluated for large clusters using QMC calculations, in 
contrast to the exact gap, which is difficult to evaluate 
directly (because it is dominated by statistical errors if 
the gap is small) . We also found that the probability dis- 
tribution of local gaps A^ (also defined using a sum rule) 
scales with the system size.^"* Defining = A^L" (where 
the exponent is determined from simulation data and is 
a « 2.8 for 5' = clusters), the distribution P(ei) is size- 
independent. Moreover, the low-energy tail of this distri- 
bution is well described by a power-law, P{ei) oc e^'^, with 
uj = 1. Analyzing the local gaps using extremal- value 
statistics, we found that the dynamic exponent should 
be related to the parameters of the local gap distribution 
according to z = a + Df/{oj + 1). Our simulation results 
satisfy this exponent relation remarkably well. The ap- 
plicability of the exponent relation supports the notion 
that the low-energy excitations involve a number oc n fi- 
nite regions (containing the effective moments), while an 
exponent a > shows that individual excitations are not 
localized (since for localized excitations the energy should 
be independent of L for large L) . The effective moments 
should be located in regions of imbalance in the number 
of spins on the two sublattices, and many moments can 
be involved in an excitation. The value of the exponent a 
reflects the way in which the weak interactions between 
the effective moments involved in a particular excitation 
decrease with increasing system size, as these moments 
become further separated from each other. 

In this paper, we report scaling results for larger clus- 
ters than previously and also compare results for clusters 
constructed in different ways. On the bipartite square 
lattice, we denote the number of sites on sublattice A and 
B as riA and n b , respectively. The ground state has spin 
S ~ \nA — nB\/2. We analyze in detail both S ~ and 
5 > clusters at the percolation point p* . We use the gap 
upper-bound A* from sum rules, as well as Lanczos exact 
diagonalization results for the excitation spectrum. For 
clusters with ground-state spin 5' > 0, we point out that 
the spectral weight entering in the sum-rule approach is 
dominated by S* — s- 5-1-1 excitations, whereas the lowest- 
energy excitations typically correspond to 5 — *■ 5 — 1. 
The quantity A* in this case describes only excitations 
where S S + 1, for which we find z « 1-5-D/ based on 
finite-size scaling. However, the lower S ^ S — 1 exci- 
tations most likely follow the same z « 2Df scaling as 
the 5 = — > 1 excitations of ua = ns clusters. We 
also discuss results for the dimer-diluted bilayer at p*, 
as well as the single layer at p < p* . For these systems, 
we observe behavior consistent with quantum rotor exci- 
tations (although other scenarios, e.g., fractonS )^^'^^ are 
also possible). 

To explain the existence of localized moments in the 
percolating cluster, we also introduce a classical dimer- 



3 



monomer aggregation model to study the purely geomet- 
rical local sublattice imbalance, which we believe is at the 
heart of this problem. The dimers correspond to nearest- 
neighbor sites that can form minimal local quantum ro- 
tors, and the monomers lead to "dangling" spins that 
are, due to local sublattice imbalance, left over after the 
maximum number of dimers has formed. The monomers, 
individual ones or groups of several of them, can lead to 
effective magnetic moments in the spin system. We find 
that the classical monomers indeed are confined within 
regions of finite size, both at and away from the perco- 
lation point. The anomalous dynamics with z k, 2Df in 
the single-layer quantum spin system at p* should there- 
fore be a consequence of localized quasi-free magnetic 
moments interacting very weakly because of the vanish- 
ing spin stiffness of the percolating cluster.^ Away from 
the percolation point, the moments can lock to the global 
Neel order of the cluster (as a single magnetic impurity 
in two dimensional is known to dc^^-^*) and do not form 
an effective independent low-energy system. 

To further investigate the nature of the excitations of 
the quantum spins and their relationship to the classi- 
cal monomers, we have also applied a projector QMC 
method in the valence bond basiai^ to directly study the 
triplet excitations of clusters with singlet ground states. 
In the valence bond basis, a triplet state can be described 
by a lone triplet bond, the location of which fluctuates 
among the background singlet bonds. We find that the 
triplet bond is indeed predominantly localized at a subset 
of the classical monomer regions. The total size of the 
excitation (i.e., the number of spins involved in it) is not 
finite, however, but grows with the cluster size according 
to a non-trivial power law. 

The outline of the rest of the paper is as follows. After 
defining the spin models and describing several compu- 
tational methods in Sec. |TT1 we present results of both 
Lanczos exact diagonalization and sum-rule QMC cal- 
culations for single-layer clusters at p = p* in Sec. IIIII 
In Sec. |IV] we discuss the distribution of spectral weight 
in the dynamic structure factor originating from excited 
states of different total spin, using Lanczos exact diag- 
onalization as well as an approximate analysis based on 
valence bond states. We discuss scaling results for per- 
colating bilayer clusters in Sec. |Vl and for single-layer 
clusters away from the percolation point in Sec. IVII The 
classical dimer-monomer aggregation model is discussed 
in Sec. IVIIl and results of the valence-bond projector 
QMC simulations of triplet excitations in Sec. IVIIII We 
conclude in Sec. IIXI with a summary and discussion. 



II. MODEL AND METHODS 

The Heisenberg Hamiltonian on a single site-diluted 
layer is given by 

H = J^5,5,S, -Sj, (J>0), (1) 

(id) 



where {i,j) denotes nearest neighbors on a 2D square 
lattice and Si = (vacancy) and Si — 1 (magnetic site) 
with probability p and 1 — p, respectively. We study 
clusters with two types of boundary conditions. In open- 
boundary L X L systems, we start with all magnetic sites 
and introduce vacancies with probability p. We study the 
largest cluster of connected magnetic sites. The num- 
ber of spins n in such clusters fluctuates and scales as 
(n) - L^f, with Df = 91/48. We also study clusters 
grown on an infinite lattice. Starting from a single mag- 
netic site, we add more sites to the cluster with prob- 
ability 1 — p by transversing along the boundary sites, 
leaving sites unfilled with probability p, but flagging each 
site as visited (so that sites assigned as vacancies are not 
visited again). This procedure terminates at random at 
some stage where all neighbors of the cluster have been 
assigned as vacancies. We only keep clusters of some de- 
sired target size n. These clusters have a characteristic 
average length (L) (defined, e.g., as their radius of gy- 
ration) such that n cx {L)^f . The two types of clusters 
will be referred to as L x L and fixed-n, respectively. In 
Ref. Ill, we only studied fixed-n clusters. Here we also 
consider the Lx L variant to check whether the finite-size 
scaling properties depend on the boundary conditions in 
the cluster construction. For p < p* , we consider only the 
L X L clusters, because the fixed-n construction rarely 
terminates at reasonably small n in this case. 

Under each type of boundary condition, we further 
consider two different ensembles of sublattice occupa- 
tions; ua = riB, in which case all clusters have ground 
state spin S' = 0, as well as arbitrary n — ha + ns (with 
the distribution give by the cluster construction), corre- 
sponding to ground state spin S = |n^ — nB|/2. The 
latter ensemble includes also the S = clusters. 

A bilayer cluster is constructed by coupling two iden- 
tical single-layer clusters with an interlayer coupling con- 
stant J±. The Hamiltonian is thus 

H = J SiSj (Sii ■ Sij + S2i ■ S2j) 

+J±^^SiSii ■ S2i, (2) 

i 

where the subscripts 1,2 refer to the two layers. Also 
in this case we can study L x L or fixed-n clusters, but, 
in contrast to the single layer, the ground state of a bi- 
layer cluster is always a singlet because each spin can be 
paired with its neighbor in the opposite layer. We con- 
sider small coupling ratios J±/J, for which the ground 
state has long-range order. ^ 

Here our main interest is in the the energy gap A 
between the ground state and the first excited state, 
which in the case of an ua — ng cluster is a singlet- 
triplet gap. For clusters with general uatTIb such that 
S — \'nA — nB|/2 > 0, the lowest excitation can have 
total spin S" = S" — 1, 5", or S" -I- 1. In addition to the 
gap, the distribution of the spin S' of the lowest-energy 
excitation is also interesting. We will also study the lo- 
calization properties of the excitations very explicitly, by 
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formulating the problem in the valence bond basis and 
carrying out unbiased quantum Monte Carlo calculations 
of S" = 1 excitations of clusters with 5 = ground states. 

To calculate the gaps, we use both direct and indirect 
(approximate, through sum-rules) estimates, using the 
methods discussed in Sees. Ill Al and ITTbI In Sec. Ill Cl we 
will introduce the valence bond QMC scheme for directly 
imaging the spatial distribution of triplet excitations. 



A. Exact diagonalization 

The most straight-forward approach is to diagonalize 
the Hamiltonian numerically in sectors of different mag- 
netization, 



(3) 



using the Lanczos method. However, for irregular clus- 
ters (without lattice symmetries to exploit for block- 
diagonalization), this can be done in practice only for 
up to n « 20 spins, due to the rapid growth of the ma- 
trix sizes with n (considering also that we have to av- 
erage over a large number — typically thousands — of ran- 
dom cluster realizations). Nevertheless, such calculations 
are very useful and give some important insights into the 
role of "dangling" spins in low-energy excitations. 

In addition to studying the level spectrum, focusing 
on a few low-lying states and calculating their total spin 
to classify the excitations, we also compute the full dy- 
namic spin structure factor (in the standard way with 
the Lanczos method, as described, e.g., in Ref . [20) : 



5(q, oj) = Y.\ (m|5^|0) \H{lo + Eo- E.„ 



(4) 



where is the Fourier transform of the spin operators; 



:^e''''-^5'|. 



(5) 



In a clean Heisenberg antiferromagnet on a bipartite lat- 
tice, the lowest excitation is a triplet at g = (7r,7r). We 
can use this wave-vector also for the diluted system, al- 
though the momentum is no longer conserved, i.e., the 
energy eigenstates |m) in @ are not classified by the 
quantum number q, but the spin operators are still 
completely well defined. We expect S{'k, tt^uj) to exhibit 
the largest spectral weight for the low-energy excitations, 
since these should involve out-of-phase fluctuations of 
neighboring spins. As we will see in Sec. IIII CI the dy- 
namic structure factor is of great utility in judging the 
validity of our sum-rule based approach for an upper- 
bound of the energy gap, which we discuss next. 



B. Quantum Monte Carlo and sum rules 

We use the stochastic series expansion (SSE) QMC 
method^! to calculate quantities which are closely related 
to the gap. An upper-bound A* to the ground state en- 
ergy gap A can be obtained using the static spin struc- 
ture factor S'(q) and susceptibility x(q) at the staggered 
wave- vector q = (tt, tt); 



A* = 2S'(7r,7r)/x(7r,7r) > A. 



This bound follows from the well-known sum-rules; 



/•oo 

/ dc^5(q,c.) = 5(q), 
Jo 

2 r^5(q,o.) = x(q), 
Jo ^ 



(6) 

(7) 
(8) 



which, in the way written here, are valid at temperature 
T = 0. In a system with a sole triplet mode (a hypothet- 
ical situation) with energy Wq, we get 2S'(q)/x(q) = Wq. 
Any spectral weight above this lowest mode will render 
the ratio larger than Wq. For a clean system, the lowest 
quantum rotor state is at q = (tt, tt) (whereas at other 
wave- vectors spin- waves are the lowest excitations). As 
we discussed above, we expect q = (tt, tt) to be the best 
choice for examining low-energy excitations also in the di- 
luted system, and we here focus exclusively on this case. 

The staggered structure factor and susceptibility can 
be efficiently calculated with the SSE method with 
" operator- loop" updates.— Using the standard defini- 
tions, for a given cluster of n sites the static staggered 
structure factor is 




^(Tr, tt) 



and the corresponding susceptibility is given by 



(9) 



x(7r,^) = -( ^(-1)*^-^- / drS., 



(r)5nO)), (10) 



where (pi — Xi + yi. Disorder averages are subsequently 
calculated for the ratio in ([6]) (where, it should be 
stressed, we first evaluate the ratio separately for each 
cluster, in order to obtain the gap bound specifically for 
each of them, and then take the average) using, typically, 
thousands of random realizations of either the largest 
cluster on L X L lattices or fixed-n clusters. 

For a clean Heisenberg antiferromagnet, A* is known^^ 
to scale with the system size as the true gap; (A*) ^ 
(A) ^ . This is because the dominant spectral weight 
is at the very lowest excitation energy — the spectral func- 
tion in the thermodynamic limit has a delta-function at 
the lowest energy, followed by a continuum at higher en- 
ergies. We expect similar spectral features in the per- 
colating cluster and suspect that A* should scale as A 
(and will show supporting numerical results in the next 
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section). At the very least, if the true power-law behav- 
ior is (A) ~ ~ n~'^/^f , then the value z extracted 
from finite-size scaling of A* must be a lower bound to 
the true dynamic exponent. Actually, in Sec. IIIII we will 
use Lanczos results for the dynamic structure factor on 
small clusters to show that the finite-size scaling of A* 
does not reflect the true lowest-energy excitations in the 
case of 5* > clusters, but all indications are that the 
sum-rule approach is valid for S* = clusters. 

We will also study an effective local (site-dependent) 
excitation gap 



1 1 

2^' 



(11) 



which is analogous to the gap bound (O but here the 
"local structure factor" is just a constant; Si = (Sf)^ = 
1/4. The local susceptibility Xi is defined as 



/3 



(12) 



Although the imaginary-time dependent correlation func- 
tion {S^{t)S^{0)) is asymptotically, for t — > oo, domi- 
nated by the lowest excitation, in practice the integral 
will be dominated by the excitation(s) which predomi- 
nantly affects the given site i. For a disordered system, 
different sites can be affected by different excitations, and 
Ai then represents a typical energy scale of excitations 
affecting spin i. 

We should note that for clusters with ground state spin 
S > 0, the grand-canonical SSE method samples over 
all magnetization sectors —S < < We there- 
fore have to subtract the static {lu = 0) contributions in 
Eqs. ([9|), (fTO]) . and p2|) arising from a non-zero ruz, i.e., 
in Eq. we subtract computed in the different 

ruz sectors and averaged over all m^. 

The SSE method operates at T > 0, but we can achieve 
the r — > limit by choosing T sufficiently low for all 
quantities of interest to converge. We use a "/? doubling" 
procedure,— in which the inverse temperature is succes- 
sively doubled until there is no longer any detectable de- 
pendence of calculated quantities on /3. Since the dy- 
namic exponent is large, the temperature T <c A has to 
be very low indeed for large clusters. As an example of 
the ultra-low temperatures required, the largest /3 we use 
for n = 512 clusters with S* = is /3 = « 5 x 10^. 
Since the simulation (CPU) time and memory usage 
scale essentially linearly in both (3 and n, these cal- 
culations are quite demanding. Fortunately, the SSE 
code for the isotropic Heisenberg model can be effectively 
parallelizedf^^ and we have run most of the simulations on 
a massively parallel computer very well suited for these 
calculations.^'* 

When studying disorder averaged static properties, the 
SSE runs for each individual cluster can be rather short. 
As long as each run is properly equilibrated (for which 
the /3-doubling procedure also helps^), the average over 
many realizations will give an unbiased estimate to any 



simple average, e.g., a spin correlation function. How- 
ever, when computing nonlinear functions involving sev- 
eral quantities, such as the ratio ([6]), the statistical er- 
rors introduce a bias. It is therefore important to collect 
sufficient statistics for the individual clusters. We have 
compared results of runs of different lengths in order to 
make sure that the results presented here do not suffer 
from significant bias effects. 



C. Valence-bond projector Monte Carlo 

To study the nature of the lowest triplet excitation of 
clusters with 5 = ground states, we apply a valence 
bond projector Monte Carlo method j^^'^^ This method 
has been described in detail in recent paper o^^i^^ and we 
here only review the elements necessary to understand 
the way we can access the triplet excitations and study 
their spatial distribution on the clusters. 

First, consider the singlet ground state |0)s, which we 
want to project out from a singlet "trial" state The 
latter has an expansion in all singlet energy eigenstates; 



I*) 



Cn\n) 



(13) 



In the standard way, if the ground state energy is the 
eigenvalue of the Hamiltonian which is the largest in 
magnitude, which can always be assured by subtracting 
a constant from H (which we assume has been done, if 
necessary), the ground state can be projected out by ap- 
plying a high power of H to the trial state; 



{-Hr\^) 

|0). 



p 



(14) 



Cl 
Co 



El 
Eo 



ID 



where we include the minus sign because normally Eq < 
0. For large P all the excited states are filtered out be- 
cause the ratios \En/Eo\ < 1. 

Valence-bond basis states are products of N/2 singlets. 



(*,j) = (Td, -i,T,)/V2, 



(15) 



where we consider the first, i, and second, j, spins to 
always be on sublattice A and B, respectively. The trial 
state is thus expressed in this over-complete basis as 



wJVy) 



(16) 



where v labels all the different tilings of the cluster into 
valence bonds, of which there are (iV/2)!, and we have 
introduced the short-hand notation \Vy) for a valence- 
bond basis state. 
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For a trial state in a ground state projector calculation, 
it is convenient to use an amplitude product state^liSS 
where the expansion coefficients are given by 



(17) 



where h{x,y) > and n^{x^y) is the number of bonds 
of size {x^y) in the configuration, i.e., the length of the 
bond is r = (a;^ +?/^)^/^. Note that it is not necessary to 
normalize the trial state. 

For the clean 2D system, the optimal amplitudes are 
translationally invariant and decay as h{r) ^ r^^f22, For 
random clusters, the optimal amplitudes are naturally 
not translationally invariant. While the average bond 
probabilities (which are related to the amplitudes) de- 
cay with r, for any given cluster there are typically some 
regions spanned by long bonds (a feature intimately con- 
nected with the low-energy physics, as we will discuss in 
Sees. IIVI and [VIIip . One could in principle optimize all 
the (X different amplitudes for each specific cluster. 
However, the effort involved in individual optimizations 
for hundreds or thousands of clusters does not necessarily 
pay off, compared to just projecting the trial state with 
a somewhat larger power P of H . We here use a very 
simple trial state with all h{x,y) = 1. 

To carry out the projection using Monte Carlo sam- 
pling, we write the S = 1/2 Heisenbcrg Hamiltonian in 
terms of singlet projection operators on all the pairs h of 
nearest- neighbor sites {i{h),j{h))\ 



Hh — H, 



and write the projection operator in (|15p as 

\b=l J r 



where 



Vr - Hh 



■ Hhr Hh 



(18) 



(19) 



(20) 



denotes the possible strings, r = 1, . . . , Nf, of the singlet 
projectors. 

When a singlet projector Hij acts on a state with a 
valence bond on the two sites the state remains un- 
changed with a matrix element of unity; we call this a 
diagonal projection. If the operator acts on a state with 
no valence bond on the two sites, then the two bonds 
(j, k) and (Z, j) connected to i,j are broken, and new sin- 
glets and (/, k) are formed. This process has matrix 
element 1/2, and we call it an off-diagonal projection. 
Thus the projection rules are; 



H,,\...i^,J)...) = \...{^,J)..■), 



(21) 



H,,|...(z,A;)...(/,j)...) = i|...(z,j)...(;,fc)...). (22) 

Acting on a component \Va) of the trial state, a string 
Vr effects a number of rearrangements (1^^ of pairs of 



valence bonds, resulting in another valence bond basis 
state which we call \Va{r)); 



Vr\Va) ^War\Va{r)). 



(23) 



Here the "projection weight" War for a combination of 
operator string Vr and state \ Va) is given by the number 
TOoff (a, r) of off-diagonal operations ([22]) in the course of 
the projection; 



Wnr = 2-"°«("''^). 



(24) 



The expectation value of an operator A can be written 

EabErl^JaWtinlVtAVrlVa) 



(A) 



EabEriy^aWb{Vt\V;Vr\Va) 
EabErl^aWhWarWuiVbimiVair)) 
EabEriy^aWbWarWtimWair)) ' 



(25) 



where Wa and Wf, are the weights computed according 
to pT]) for the bonds in the states 1^4) and (Vf,| in the 
expansion (|16p of the trial ket and bra ^('I'l states. 

The sampling weight to be used in Monte Carlo calcu- 
lations of (pS)) is 



W{a,b,r,l) ^ WaWbWarWu{Vb{l)\Va{r)), 



(26) 



where the overlap of the two projected states is given by 

mi)\Va{r)) = 2^-^/2, (27) 

where A^o is the number of loops formed when the bond 
configurations of the states |Va(r)) and (Vb(Z)| are su- 
perimposed (forming the transposition graph^^). Simple 
sampling procedures for the operator strings and trial 
state bonds are described in Refs. [T9II26I . More efficient 
sampling methods have been developed recently^^^ which 
we use but do not discuss here. 

For the purpose of the present paper, the most inter- 
esting aspect of the valence bond projector scheme is the 
fact that we can easily extend the scheme to also study 
a triplet state. A trial wave function in the triplet sec- 
tor can be expressed in the overcomplete basis of a lone 
triplet bond among A^/2 — 1 singlets. We denote a zero- 
magnetization triplet by square brackets; 



and expand the triplet trial state as 



N/2 



N/2 



(28) 



N/2' j N/2)) 



(29) 



where the normalization is again irrelevant. Here we use 
the same expansion coefficients — the amplitude products 
(fT7)) — as in the singlet trial state. Note that for a clean 
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system, the singlet state ^TE\\ has momentum k = (0,0), 
whereas the triplet has k = (tt, tt). These are the 
known momenta of the lowest states in the two spin sec- 
tors (with the triplet being the lowest member of An- 
derson's tower of quantum rotor states^^). The wave- 
function signs corresponding to (fT6|) and ([29|) should be 
correct for the lowest singlet and triplet states also for 
a diluted system, since all conditions for Marshall's sign 
rule [which corresponds to all positive expansion coeffi- 
cients in Eq remain valid. 

When acting on a triplet bond, the singlet projector 
Hij destroys the state, while the action between a singlet 
and a triplet bond are very similar to the pure singlet 
rules P^ . The two triplet rules are 

H.j\...[hj]-)=0, (30) 
7J.,|...[z,fc]...(;,j)...) = i|...(z,j)...[Z,fc]...). (31) 

In the projector method, it is straight forward to convert 
one bond of the singlet trial wavefunction into a triplet 
and trace its evolution. The triplet states that survive 
after all P operations [i.e., that are not destroyed by a 
diagonal operation ((50)) ] are used to measure properties in 
the triplet sector. To measure triplet expectation values, 
we have to project triplets like this both in the bra and 
ket in the triplet version of (|25p. We also have to check 
the overlap ([26]) of the surviving triplet states. One can 
show that the two triplet bonds have to be in the same 
transposition-graph loop in order for the overlap to be 
non-zero, and it is then equal to the singlet overlap (P7)) . 
For surviving pairs of triplets, the weight of the triplet 
configuration is the same as that of the original singlet 
one. One can therefore sample the configurations in the 
singlet sector, and carry out measurements with all the 
surviving triplets without reweighting. This is one of the 
strengths of the valence bond projector method. 

There can still be problems with this approach, be- 
cause the number of surviving triplets decreases with the 
projection power P [because the probability of a triplet 
to be destroyed by a diagonal triplet operation ((30)) in- 
creases). It helps considerably that the starting trial 
state can have the triplet at N/2 different locations, in 
both the bra and the ket state, and as long as one pair 
out of the total of combinations survives (and 

gives non-zero overlap), we can collect statistics. One 
can carry out the summation over triplet locations m in 
((^5)) efficiently, without introducing any additional factor 
N/2 in the computational effort, in a single traversal of 
the operator sequence. 

In some cases it can still happen that the triplet quan- 
tities of interest have not converged well to the P — > oo 
limit before the triplet survival probability becomes too 
low to be useful. This is not a serious problem in the 
present application, although an extrapolation to infinite 
P based on several calculations with reasonable triplet 
survival probability is necessary to ensure that the results 
represent the lowest triplet. An exponential asymptotic 
convergence can be expected based on Eq. (fT5)) . 




FIG. 1: (Color online) Results for two different clusters, visu- 
alized with color scales, of the QMC sum rule approximation 
of local gaps Ai (left) and valence bond projection calculation 
of the triplet density pi (center). To the right, the clusters are 
shown covered with dimers (pairs of spins enclosed by ovals) 
and left-over monomers (black circles). The black circles in- 
side ovals indicate other possible locations of the monomers, 
corresponding to alternative maximal dimer coverings (which 
here always corresponds to two left-over monomers). 



We will discuss the spatial distribution of the triplets. 
The surviving triplet states have the triplet bond located 
at two particular sites (which can be different in the ket 
and the bra, and we do the measurements in both of 
these states). In a random system, the average triplet 
density will not be uniform and provides a very concrete 
measure of the localization properties of the lowest triplet 
excitation. 

Note that the distribution of the to^ = triplet bond 
is equivalent to the magnetization distribution in a state 
with ruz — 1, which could also be studied using the SSE 
method at low temperatures (e.g., by including a weak 
magnetic field^). However, the valence bond states also 
contain other relevant information, e.g., the statistics of 
the length of the triplet bond, which can only be accessed 
in the valence bond basis and which will be useful for 
analyzing the nature of the excitations (as we will do in 
Sec. (Villi)- 



D. Examples 

Having introduced the technical aspects of all the 
methods, we now present illustrative results for two small 
clusters. This will help to clarify the subsequent analysis 
and discussion of results for larger clusters. 

The local gaps and the triplet density pi are visu- 
alized for two different clusters in Fig. [TJ Here the color 
scales were created separately for the two clusters, with 
the minimum and maximum values for each quantity on 
a particular cluster corresponding to the extrema of the 
scales shown (and, thus, the plots should only be used to 
examine the variations within the clusters, not compar- 
ing the values for the two clusters). 

The two clusters differ qualitatively in a way which 
is directly related to our arguments pertaining to a low 
energy scale. The lower cluster can be completely sub- 
divided into pairs of nearest-neighbor sites (dimers, rep- 
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resented by ovals), whereas the upper one has two "dan- 
ghng spins" left (monomers, shown as black circles out- 
side ovals) after the sites have been paired up as much 
as possible. The pairing into dimers is not unique — the 
black circles inside ovals show all other possible monomer 
locations for this cluster. In all cases there are two 
monomers in two separate regions. The classical dimer- 
monomer aggregation model discussed in Sec. IVIII con- 
tains the statistics of the distribution of the monomers. 
Our main argument is that the presence of monomer sites 
leads to small gaps, i.e., a large dynamic exponent. For 
the two clusters shown, the exact gaps are 0.039J and 
0.276J, respectively, for the cluster with and without 
monomer sites. The gap upper bounds A* are 0.076J 
and 0.35J. While in particular the former is quite far 
from the exact result (in a relative measure), the differ- 
ence between the two clusters is still large. 

Large clusters are likely to have dangling spins, and 
the top cluster in Fig. [T] is therefore the more interesting 
case. One can clearly see a strong correspondence be- 
tween small local gaps and large triplet density, and they 
both coincide very well with sites where monomers can 
be located. Although this is in accord with the notion 
of monomers leading to finite regions of spins affected by 
the excitation, these clusters are clearly too small to give 
any meaningful quantitative insights into the localization 
properties of the triplet. 

In the following four sections we will carry out quan- 
titative scaling analyses of the gaps in different types of 
clusters, while further discussion of the monomer and 
triplet distributions will be postponed to Sees. IVIII and 
IVIIIl respectively. 



III. SINGLE-LAYER GAP SCALING AT THE 
PERCOLATION POINT 

We here discuss the distribution of exact gaps obtained 
with the Lanczos method, as well as SSE QMC results for 
the gap upper bound and local gaps. First we consider 
5 = clusters {ua = ns), and then arbitrary S. 



A. Clusters with singlet ground state 

Fig. mja) shows the probability distribution of the log- 
arithm of the exact gap A of n = 16 clusters obtained 
using 4x 10* samples. We also show results for the upper- 
bound A* for clusters of the same size, obtained from 
SSE calculations for 6 x 10^ different clusters. We pre- 
sented these results in Ref. [13 and here re-graph them in 
a different way for added clarity. The A* curve is visibly 
shifted up in energy relative to the A distribution (with 
the average A*/A « 1.5), but the shapes of the two 
curves are remarkably similar. The two-peak structure 
is related to the "danghng spins" discussed in Sec. IIIDI 
The large-gap peak originates almost exclusively from 
clusters that can be completely partitioned into nearest- 




FIG. 2: (Color online) Distribution of the singlet-triplet gap 
A and its upper bound QMC estimate A* for n — 16 clusters 
(a), and the A* distribution for n = 32 (b) and n = 64 (c). 



neighbor dimers, whereas the low-gap peak corresponds 
to clusters with dangling spins (monomers). Clearly, as 
the cluster size grows, it will be less and less likely to find 
clusters with no monomers, and the weight of the high- 
energy peak should therefore gradually diminish and be 
absent for large clusters. The relative size of the large-A* 
peak is indeed much smaller in the n = 32 distribution 
graphed in Fig. [IJb). In the L = 64 histogram, shown 
in panel (c), only a single peak can be discerned (with 
only a weak tail suggesting some remaining contributions 
from no-monomer clusters). 

Fig. El shows the size dependence of the disorder aver- 
aged (A*) on log- log scales for both fixed-n (top panel) 
and LxL (bottom panel) clusters. We also show the typ- 
ical values (A*)t, obtained by averaging ln(A*) for the 
individual clusters. While the typical and average values 
do not exactly coincide, for large systems they scale in 
the same way. Linear fits to the (A*)t data on the log- 
log scales gives z = 3.6 ± 0.1 for both types of clusters. 
Here the estimated error reflects the purely statistical er- 
rors of the line fits in combination with small variations 
depending on what range of system sizes are included. 

As shown in Figs. HUa) and[5fa), for fixed-n and L x 
L clusters, respectively, not only do the averages and 
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FIG. 3: (Color online) Finite-size scaling of the average (A*) 
and typical (A*)t gap upper-bound for S = (ua = ris) 
clusters. The top and bottom panels show results for fixed-n 
and L X L clusters, respectively. The lines correspond to the 
scaling expected with dynamic exponent z — 3.6 (i.e., the size 
dependence is ~ n~^^^^ and ~ , respectively, for the two 
types of clusters). 
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FIG. 4: (Color online) Distribution of the logarithm of the 
scaled gap upper-bound e = A'n'^^^f (a) and local gap bound 
ei = Am'^^^f (b) for fixed-n clusters with ground state spin 
S = 0. The exponents are indicated in the panels. The solid 
lines correspond to small-gap exponent ui = 1 and the curve 
in (a) is a Frechet form. 



typical values of A* scale with the system size, but the 
entire distribution can be collapsed onto a common size- 
independent curve, by scaling the gap estimates with the 
cluster size. We define the scaled gap upper-bounds for 
the two types of clusters according to 

_ J A*L^, (fixed-n clusters), , , 

^ ^ \ /\*n'/°i, {LxL clusters). ^'^ ' 

As can be seen in the figures, the small-gap side of the 
distribution of ln(e) is very well described by a power law; 
P[ln(e)] cx 6^^+^ with cj = 1. This distribution of the 
logarithm of e corresponds to a probability distribution 
P(e) ~ e'^ for the scaled gap e itself [since the differential 
(iln(e) — de/e]. 

We next consider the local gap estimate A^, i.e., the 
inverse local susceptibility p2|) . Measuring this quantity 
at each site, we define size-scaled local gaps; 

_ J AiL^, (fixed-n clusters), , . 

^'~\A,n"/^/, (L X L clusters). ^ ' 

The probability distributions of In(ei) for different clus- 
ter sizes, based on several hundred clusters of each size, 
collapse onto each other for a suitably chosen a « 2.8, as 
shown in Figs. IDJb) and[n)Jb) for the two types of clus- 
ters. The small-gap tails of the distributions are again 



very well described by a power law; P{ei) ^ ef , with the 
same a; = 1 as for the scaled "global" gap bound A*. 



B. Extremal- value analysis 

In Ref . [3 we used extremal value statistics'^^ (in a way 
generalizing a treatment of localized excitations by Lin 
et ali^) and found a relationship between the exponents 
z, a, and lo. For completeness, we repeat and further 
clarify our arguments here. 

Our hypothesis is that, for a large cluster of size n, 
there is a number oc n of regions of sublattice imbal- 
ance. These regions act as localized magnetic moments, 
which interact weakly with each other through the mag- 
netically inert parts of the percolating cluster. The ex- 
citations of this effective low-energy system of coupled 
moments are not localized because several distant mo- 
ments can be involved. It is then natural to expect some 
size dependence of the local gaps, due to the dependence 
of the effective interactions on the distance between the 
moments involved in a low-energy excitation, combined 
with the increasing distance (on average) between these 
moments with increasing cluster size. We posit that this 
size dependence can be captured by the single exponent 
a in Eq. 
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FIG. 5: (Color online) Probability distribution of the loga- 
rithm of the scaled gap upper-bound e = A*L^ (a) and local 
gap bound — A^L" (b) for Lx L clusters with ground state 
spin 5 = 0. The solid lines correspond to lj = 1 and the curve 
in (a) is a Frechet form. 



The actual finite-size gap A for a given cluster should 
correspond to the smallest of the local gaps Ai for that 
cluster, for which we use the notation Amin. Of course, 
the local gaps that we measure are only approximations; 
one cannot unambiguously define a local gap in an inter- 
acting system. Nevertheless, A^ reflects the local distri- 
bution of spectral weight, and there should be some site 
i within the regions affected by the lowest excitation for 
which Ai = Ajjiin « A (and Amin > A). In our numer- 
ical analysis, A is approximated by the bound A*, and 
we expect A^in ~ A*. Examining the numerical data, 
we indeed find a very strong correlation between the two 
quantities, as shown in Fig. [5] for fixed-n clusters. Here 
it can be seen that Amin is typically 1.5 — 2 times larger 
than A*, which reflects larger spectral weight above the 
true lowest excitation energy in the local dynamic struc- 
ture factor Si{oj) than in S'(7r,7r, w). It should be noted 
that Amin < A* is allowed within the sum-rule approach, 
although Amin > A has to hold strictly. 

We now assume that there is a number M cx n of differ- 
ent local scaled gaps and investigate the consequences 
of this in light of the scaling behavior found above. We 
assume a probability distribution P{ei) — Ae^ for some 
window of small e (where ^ is a constant and we consider 
a more general case than just a; = 1 extracted from the 
finite-size scaling of the data). We derive the probability 
distribution for the smallest scaled local gap PAf(emin) 
for large M ~ L^s using extremal-value statistics. 
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FIG. 6: (Color online) Correlations between the smallest local 
gap Amin and the gap upper bound A* for clusters of size 
n = 64 and 128. Each data point corresponds to Monte Carlo 
results for a randomly generated ua ~ ?^s cluster. The line 
shows the ideal (single-mode) case of complete equivalence of 
the two estimates of the finite-size gap. 



We should clarify why we assume M oc n for the num- 
ber of local gaps, instead of just M = n, which is the 
actual number of different numerical values A^ that we 
compute for a given cluster. The distinction will not 
matter in the analysis, but it has an important physical 
significance. In our scenario, a cluster consists of regions 
with localized moments, which participate in the low- 
energy excitations, as well as inert parts which have only 
high-energy excitations. The form of the probability dis- 
tribution P{ei) = At^ should only hold for sites i within 
the moment regions. It is then important in our analysis 
that also the number of such sites scales as n (although 
one could also generalize to M ~ n''' with 7 < 1, but 
the consistency of our analysis with M cx n will show 
that this is not necessary). Later, we will provide more 
concrete proof that the moment regions arc finite and 
the total number of spins belonging to moments grows 
linearly with n. 

Note also that in reality the distribution of local gaps 
must be cut off (equals zero exactly) below some very 
small value for a given finite cluster size. However, this 
should not affect the results of the analysis to follow, 
because also the assumed power-law probability is very 
small below such a threshold. We thus expect the results 
derived below to be valid within some significant window 
of scaled gaps e. 

We denote the probability of finding a local gap at an 
arbitrary chosen site (within one of the moment regions) 
smaller than some value x by P^{x). It is given by 



P<(x) 



A 



(34) 



If one of the scaled gaps ej is the smallest and has the 
value e, then all the other (M— 1 different) values ei, i j 
must be larger than e. The probability of these M — 1 
values being smaller than e is [1 — P<(e)]*^~^. Since any 
of the M values could be the smallest one, we get a factor 
of M, and finally the distribution of the ej value is given 
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by P{€j)- Thus, the distribution of the smallest scaled 
local gap is 

PM{e) = MP{e)[l~P^{e)r-\ (35) 
which for small e also can be expressed as; 

PMie) = - P<ie)f - "l^^"'''^^^^- (36) 

Using Eq. ([M]) here gives the Frechet distribution;^ 

Pf{u) = Au'^expi^Aiu; + (37) 

where u = uoe. Thus, the probability distribution of the 
scaled global gap should be governed by the same expo- 
nent Lu as the scaled local gaps. The Frechet form can 
indeed be fitted to the A* data in Figs. Hl^a) and[5Ka), 
with the same exponent w = 1 as in the local-gap (b) 
panels, but only in the small-gap region. One cannot ex- 
pect the Frechet distribution to work for large gaps, since 
the local gap distribution we started from is linear only 
in the small region (and, as discussed above, we expect 
the large-gap part of the distribution to be dominated by 
excitations of the magnetically inert cluster regions with- 
out moments). The fitted forms in Figs. [DJa) and [Hl^a) 
are therefore also not normalized. Nevertheless, it is en- 
couraging that the data is in agreement with the result 
that both the local and global gaps should scale with the 
same exponent, which here is w = 1. 

Let us now use the distribution p4|) in a different 
way. Since we assume that there are M oc n local gaps, 
the typical smallest gap should correspond to P< (x) for 
which X = M~^, i.e., x cx L~^f . This gives Amin oc 
j^-a-Df/{uj+i)_ g-j^(,g ^^.^ should equal A, and, by defi- 
nition, A cx we arrive at the following relationship 
between the three exponents; 

z = a+-^. (38) 

This generalizes the relation z = Df/{uj + 1) used as a 
criterion for a localized excitation by Lin et al.'^^ to ex- 
citations originating from two or more finite entangled 
regions distributed over the cluster. With our numeri- 
cal values from the finite-size scaling above, z « 3.6 and 
w = 1 (the latter of which is not based on a fit, but is a 
value consistent with all our results), we obtain a ~ 2.65, 
in very reasonable agreement with the value a w 2.8 ob- 
tained in Figs. HUb) and ^h) from the scaling of the 
data for the fixed-n and L x L clusters. The applicabil- 
ity of the exponent relation (|38p provides strong support 
to our hypothesis of "globally entangled local moment 
excitations" . 

C. SSE results for general-S clusters 

We now turn to clusters with no restriction on the sub- 
lattice occupations and ng in the generated ensem- 
ble. Scaling results for the global and local gaps of Lx L 
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FIG. 7: (Color online) Distribution of the scaled gap upper- 
bound (a) and local gaps (b) of Lx L clusters with no restric- 
tion on the sublattice occupation numbers ua and ub- The 
solid line shows the asymptotic small-gap behavior expected 
with uj = 1. 



clusters are shown in Fig. [71 The finite-size scaling of 
the average and typical values of A* are shown for both 
fixed-n and L x L clusters in Fig. [51 We obtain a « 2.1 
and z « 2.8 for both cluster types. These exponents 
differ significantly from the ones obtained previously for 
the ensemble including S* = clusters only. In particular, 
z « 1.5-D/, whereas the S = clusters gave z « 2Df. 
The exponent relationship (|38p still holds approximately, 
albeit with somewhat larger deviations than in the S = 
case. The small-gap behavior remains consistent with the 
exponent w = 1 in all cases. 

We believe that the much smaller exponents z, a are 
due to a failure of the sum rule approach to capture the 
true low-energy states for 5* > 0. To demonstrate this, 
we next investigate the dynamic structure factor ([?]). 



IV. SPECTRAL WEIGHT DISTRIBUTION 

We will investigate how the spectral weight of the dy- 
namic structure factor is distributed among different spin 
sectors of the excited states in (U). Acting on the ground 
state with the q = (tTjTt) spin operator (O, or the cor- 
responding X 01 y components, on one of the (25* -I- 1) 
degenerate ground states of spin S results in states with 
spin S and 5 ± 1 . This well known selection rule for the 
dynamic structure factor ^ can be easily demonstrated 
in the valence bond basis. Here we do this as a prelude to 
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FIG. 8: (Color online) Finite-size scaling of the average and 
typical gap upper-bound A* for the ensemble with unre- 
stricted ua and ub ■ The upper and lower panels show results 
for fixed-71 and Lx L clusters, respectively. Lines correspond- 
ing to a dynamic exponent z = 2.8 are shown with all the 
data sets. 



discussing the distribution of the spectral weight among 
the three sectors of final spin for clusters with ground 
state S > 0. 



A. Selection rules 

We consider an extended valence bond basis with an 
arbitrary number of = triplet bonds (|29p in addition 
to singlet bonds, for a state with total niz = (hence the 
number of spins, N, is even). Later, we will consider also 
rUz 0. The standard valence bond basis for ua — ns 
is restricted to bipartite bonds onlyi^l Here, for ha ^ ub 
and rriz — 0, we require a maximal number of bipartite 
bonds, i.e., if the total sublattice imbalance is defined as 
Aab = — "-bI/S, there will be ni, — N/2 — Aab bipar- 
tite bonds and Uc = Aab bonds connecting sites on the 
same sublattice (with all such pairs either on the A or B 
sublattice, depending on which sublattice has the larger 
number of sites). This basis is clearly overcomplete. A 
state with two triplet bonds and one non-bipartite bond 
is illustrated in Fig. [n{a). 

First, let us discuss the relationship between the total 
spin S and the number of triplet bonds. A state with 
rit triplets does not have fixed spin when rif > 1 (while 
for rit = and 1, the state has fixed 5 = and 1, re- 
spectively). According to the rules for addition of angular 




FIG. 9: (Color online) Valence bond states on a cluster with 
sublattice imbalance Aab = 1 (requiring one non-bipartite 
bond — here the top one). Open and solid circles indicate the 
two sublattices. Singlet and triplet bonds are shown as solid 
and dashed lines, respectively, (a) shows a state with two 
triplet bonds and rriz — 0. In (b), there are two unpaired up 
spins and lUz — 1. 



momenta, one might at first sight suspect that nt triplets 
could be used to form states with S ^ 0,1, ... ,nt. How- 
ever, consider the operator Z which inverts all the spins; 



z\s!,s^ 



• ^n) — I ^1: 



-s 



(39) 



which is a special case of a rotation in spin space. Since 
the total magnetization ruz = 0, a state with fixed S is 
also an eigenstate of this operator, with eigenvalue z — 
±1. Since a triplet pair (bond) is even under Z while 
a singlet pair is odd, the eigenvalue z of a state with 
a fixed number rit of triplet bonds is z = (—1)^/^^"'. 
Thus, in order to construct a state with fixed S (fixed z), 
one cannot mix valence bond states with even and odd 
number of triplets. Since the minimum number of triplets 
required to construct a state with fixed spin is nt = S, 
we conclude that the triplet numbers that can be mixed 
are nt & {S, S + 2, . . . N/2}. This, in turn, implies that a 
state with fixed number of triplets is a linear combination 
of states with S & {0/1, ... ,nt — 2,nt}, where the lower 
limit or 1 applies for even and odd S, respectively. 

Next, we let the q = (jr, tt) spin operator ^ act on a 
given valence bond state with nt triplets. We can write 
the operator in a way tailored specifically for the state 
under consideration; 



5f 



1 



.b=l 



c=l 



(40) 

Here the subscripts i{b) and j{b) refer to two sites con- 
nected by a bipartite valence bond b, and k{c),l{c) de- 
notes a pair of sites on the same sublattice, connected 
by a non-pipartite bond c. The bonds can be singlets or 
triplets, and the possible outcomes when operating with 
one of the terms are; 



[ib.ib]--) = \...{ib,ib)- 
{icjc)---) = 0, 
(%)+%))l-[*c,Jc]...) =0. 



i^i{b) 



S]{b))\- 



(41) 



Thus, operating with the full S^^, we obtain a linear 
combination of states with nt + 1 and nt — 1 triplets. 
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Extending this result to the case of a fixed-S' state l^'s), 
which is a Unear combination of states with different rit 
(all even or all odd), we can think of the triplet bond cre- 
ated or destroyed in each term [with the operator ([^0]) 
written in the appropriate way for operation on each 
term] as adding or subtracting a spin 1 to or from a 
spin S. Then, considering also that even and odd S cor- 
responds to mixtures of even and odd n^, respectively, 
we conclude that the state S^^^\^s) is a mixture of only 
5" ± 1 states (which is also consistent with the fact that 
for S* = ground states, the spectral weight is exclusively 
due io S = 1 excitations). 

In order to respect the spin-rotational invariance when 
using the z-component operator ^ in the dynamic 
structure factor for S* > 0, we also have to consider 
non-zero ruz- Some of the spins are then not paired up 
into valence bonds. In a minimal basis mixing valence 
bonds and spins, there are 2mz unpaired up or down 
spins for niz > and iriz < 0, respectively. The un- 
paired spins cannot be restricted to the same sublattice, 
so now the basis consists of the unpaired spins at arbi- 
trary locations, a maximal number of bipartite bonds on 
the remaining locations, and the rest of the sites covered 
by non-bipartite bonds. An example of such a state is 
illustrated in Fig. [H^a). 

In Eq. (j40p rif, and ric are the number of bipartite and 
non-bipartite bonds in a given basis state and rif, -|- n„ = 
n — rriz ■ We now also have to add a sum over the 2toz 
unpaired spins. It is then clear that ttI^s) will contain 
also a spin-S* component, arising from this added sum, in 
addition to the S" ± 1 components (which can be argued 
for in analogy with the = case). Some, but not all, 
of the corresponding spectral weight in the spin S sector 
is at w = 0, as discussed in Sec. IIIBl Averaging over all 
TTiz — —S,...,S— 1, S, it is also clear that the amount 
of S ^ S spectral weight should increase with S, as it is 
zero for S* = and the relative weight of the operations 
on unpaired spins increases with to^. 



B. Results for small clusters 

We now turn to numerical results for the dynamic 
structure factor. Investigating small clusters with the 
Lanczos method, we have found that the lowest excita- 
tion of a cluster with ground state spin S almost always 
has spin S — 1, whereas the dominant spectral weight 
arises from a state with 5-1-1. An example of this be- 
havior is shown in Fig. Iior a) for a cluster with ground 
state spin S — 3. The spectrum is dominated by a large 
contribution from an 5* = 4 state at cj/J w 1. However, 
there are numerous very small contributions from S — 2 
states below this peak, including the lowest excitation at 
Lu/J fa 0.05. In this case the sum rules give a bound 
A* very close to the energy of the lowest 5 = 4 state, 
and, thus, differs from the true gap A by a factor of 20. 
In contrast. Fig. [TUf b) shows results for a cluster with 
S = ground state. Here there are of course no exci- 
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FIG. 10: (Color online) Dynamic structure factors of two 20- 
site clusters with ground state spin S — 3 (a.) and S = 
(b). The cluster shapes are drawn in the panels. The delta- 
functions in Eq. Q are represented by vertical lines of length 
equaling the spectral weight. The symbols on top of the line 
indicate the spin of the corresponding excited states relative 
to the ground state S. In (b), there is no S — > S spectral 
weight; the circles only indicate the locations of such states. 



tations with 5 — 1, and all the spectral weight is in the 
5-f 1 = 1 channel. Moreover, the dominant weight origi- 
nates from the lowest excitation. The sum rule approach 
here gives a bound reasonably close to the true gap. 

We further examine the statistics of the gaps corre- 
sponding to excitations with 5 ± 1 for n = 20 clusters 
with ground state S. In Fig. [TT] we show histograms 
based on several hundred clusters with 5=1,2, 3, along 
with results for 5 = clusters for comparison. We can 
see that the distribution of the 5-1-1 excitations is peaked 
at higher energies than the 5—1 ones, and the distance 
between the two distributions grows with 5. As we dis- 
cussed in Sec. IIII Al the distribution of singlet-triplet ex- 
citation gaps is double-peaked for small systems, with 
the upper peak diminishing as a function of the cluster 
size. In Fig. [11] the lower part of the 5 = 0^1 distri- 
bution is located below the 5—1 distributions for 5 > 
clusters. It appears plausible from these results that the 
5^5 — 1 and 5 — > 5 -I- 1 gaps can have different scaling 
properties. 

It is also clear from these calculations that the sum 
rule approach for 5 > clusters does not reflect the 
true smallest gaps, which are due to 5 — 1 excitations, 
but instead reflect the distribution of spectral weight of 
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FIG. 11: (Color online) Distribution of the energies of the 
lowest excitation with S = S±linn = 20 clusters with 
ground state spin S. Results for S — 1,2,3 are shown in 
panels (a),(b),(c). In all panels, results for the lowest triplet 
excitation of S = clusters is shown for comparison. 



S' + 1 excitations. The quantity A* therefore has a dif- 
ferent meaning, which can still be physically relevant be- 
cause many experimental techniques probe S{q,uj) di- 
rectly, e.g., neutron scattering and nuclear magnetic res- 
onance. These experiments should observe low-energy 
dynamics corresponding to z w 1.513/ , according to our 
results in the previous section. The most plausible sce- 
nario is that the lowest S —I excitation energies, for large 
clusters and typical ground state spin (which is of the or- 
der ^/n), scale with the same dynamic exponent z « 2_D/ 
as the triplet excitations of 5 = clusters (which we will 
argue further also in the next section). While their low 
spectral weights would make them difficult to observe 
in measurements sensitive to S{q,uj), they are of course 
still relevant for thermodynamic properties such as the 
specific heat. 




FIG. 12: (Color online) Valence bond states corresponding 
closely to true eigenstates of a 6-site cluster with ground state 
spin 5=1. Solid and dashed bonds correspond to singlets 
and triplets, respectively, (a) is the ground state, (b) the 
lowest S — excitation, and (c) is obtained from (a) by acting 
on it with (which is a good approximation to the lowest 
S — 2 excitation). 



C. Valence bond theory 

We now address the important issue of why the S 
5 — 1 contribution to the spectral weight is so small. We 
will argue that this is, in fact, consistent with our scenario 
of the low-energy excitations being due to effectively iso- 
lated magnetic moments. To illustrate this point. Fig. [12] 
shows valence bond states for a 6-sitc cluster with ground 
state spin S = 1. The state in (a) is constructed as 
an approximate ground state based on the notion that 
triplet bonds should be predominantly located in regions 
of sublattice imbalance. This cluster has two "dangling" 
spins, which we take at maximum separation. For the 
two singlet bonds, we construct a symmetric combina- 
tion (which corresponds to the true ground state of the 
Heisenberg model on the four sites in isolation) . It is now 
natural to assume that the lowest excitation corresponds 
to converting the triplet bond into a singlet, as shown in 

(b) . Diagonalizing the hamiltonian exactly, we find that 
these simple states indeed are good approximations to the 
eigenstates; the overlap of (a) with the true ground state 
is 0.814, while the overlap of (b) with the lowest S — 
state (which is the lowest excited state) is even larger, at 
0.973. On the other hand, if we act with S^^^ on state 
(a) , as explained above with the spin operator written in 
the form (|40|) . we obtain the state shown in Fig. [l^c). 
This state mixes S — and S — 2 states, and its overlap 
with the actual lowest S = 2 state is 0.789. The overlap 
of (c) with the approximate S = state (b) is exactly 0, 
and the overlap with the exact lowest singlet also van- 
ishes. In the case of the (b),(c) overlap, it is immediately 
clear that it is zero because of their different states of the 
long bond. The states also differ in the quantum number 
related to a 180° rotation of the cluster; (b) is odd and 

(c) even under this symmetry transformation. The true 
ground state is also odd, which explains why the over- 
lap with state (c) is exactly 0. This latter property is 
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of course particular to this symmetric 6-site cluster. In 
general, for a less symmetric larger cluster with two dan- 
gling spins, we would expect some small overlap between 
Tries') ^^"i ttie lowest — 1 state, because the triplet 
will not be exactly localized at only two sites 

Based on the above example, we can understand that, 
in general, S* > ground states contain some triplet 
bonds connecting non-bipartite sites. The lowest exci- 
tation should normally have spin S —\ and closely corre- 
spond to converting one triplet bond into a singlet. On 
the other hand, acting with the spin operator one ob- 
tains a linear combination of 5 — 1 and S* -t- 1 states with 
an additional triplet bond, and the overlap of the S — \ 
component with the low-energy states with this spin is 
low (because of the differing singlet/triplet state of one 
non-bipartite bond). The lowest S — 1 excitation should 
thus be very similar to the excitations we have argued for 
in the case of the S — \ excitations of a singlet ground 
state, which essentially corresponds to promoting a long 
singlet (between two moments, which can be located far 
away from each other) into a triplet. For an 5 > clus- 
ter we instead demote a long triplet bond into a singlet. 
This similarity also suggests that the true dynamic ex- 
ponent (giving the scaling of the lowest energy, not the 
dominant spectral weight) in the case of S" > clusters 
should be the same z k, 2Df that we have found for the 
S = clusters. 
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FIG. 13: (Color onine) Properties of the gap upper-bound 
A* for dimer-diluted bilayer systems at inter-layer coupling 
g = 0.01 at the percolation point, (a) shows the scaling of the 
full probability distribution with z = 1.7, while (b) shows the 
size dependence of the average and typical values, along with 
lines corresponding to the asymptotic behavior with z — 1.7. 



V. BILAYER MODEL AT p" 

Our hypothesis for the low-energy excitations is that 
they are due to effectively unpairable spins on the per- 
colating cluster. To test this hypothesis further, we con- 
sider a case where there are no such spins; the bilayer 
Heisenberg antiferromagnet with "dimer dilution", i.e., 
two identical clusters coupled through a nearest-neighbor 
inter-layer coupling Jj_ = gj. The hamiltonian for this 
system was already written down in Eq. ((2). Its static 
properties were studied in Refs. [llllllll. The percolating 
cluster remains ordered at T = when the coupling ratio 
g < 0.1, whereas for larger inter-layer couplings the clus- 
ter is quantum disordered. Here we consider g = 0.01; 
well inside the ordered regime. One might then expect 
the quantum rotor picture to be valid, as has been ar- 
gued also based on field theoretical considerations,^ and, 
thus, the dynamic exponent should be z = Df « 1.89. 
Scenarios, involving "fractons" are also possible 

Fig. 1131 shows scaling results of the kind we previously 
discussed for the single layer. The peak of the proba- 
bility distribution of the bound A* for different cluster 
sizes L (the largest cluster of diluted Lx L lattices) coin- 
cides when scaled with and z « 1.7. This exponent is 
slightly smaller than Df, but considering statistical un- 
certainties of several percent and effects of subleading size 
corrections, z = Df is plausible, in contrast to z ~ 2D/ 
in the single layer. Note that the data in Fig. [T^ a) do 
not collapse onto a single curve as clearly as in the single- 



layer plots m and [51 The scaled gap distributions instead 
appear to become narrower with increasing L. This may 
be due to self-averaging following from the global nature 
of quantum rotor excitations. The local gap distribution 
(not shown here) also does not scale well with L. 



VI. SINGLE LAYER AWAY FROM THE 
PERCOLATION POINT 

An interesting question is whether the small energy 
scale of the single-layer clusters at p* survives also away 
from the percolation point. We here examine L x L 
systems diluted at p < p* , again studying the largest 
cluster for each dilution realization (which now is two- 
dimensional; (n) L^). We only consider clusters with 
ground state spin 5 = 0. 

Fig. [T3] shows results for the gap upper-bound at 
p = 0.3. For the largest few sizes the data are consis- 
tent with power-law scaling corresponding to z = D ~ 2 
(with a statistical error of sa 10%); very different from the 
behavior at p* . Given our scenario for the excitations ex- 
actly at p* , the much smaller z away from p* is either an 
indication of the moment regions not existing, or their 
mutual effective couplings (or their couplings to the rest 
of the cluster) being much stronger, thereby invalidating 
the picture of an effective low-energy subsystem. We still 
expect regions of sublattice imbalance away from p* , as 
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FIG. 14: (Color online) Scaling properties of A* of single- 
layer clusters at dilution fraction p = 0.3. For L > 28, the 
distribution can be collapsed with dynamic exponent z = 1.8 
as shown in (a). For smaller sizes, there is a cross-over behav- 
ior, as shown in (b) for the average and typical A*. The two 
lines correspond to 2 = 2 (a possible asymptotic value) and 
z' — 2.41 (a pseudo-scaling exponent in a cross-over regime). 



we will discuss further in the next section. It may not be 
surprising, however, that the moments associated with 
these are not weakly coupled, because for any p < p* 
the largest cluster has a finite spin stiffness (also in the 
thermodynamic limit), whereas exactly at p* the stiff- 
ness vanishes (although the cluster is still ordered) The 
order is thus much more robust, and as a consequence 
all the effective moments at p < p* may be locked to 
the global Neel vector and cannot be regarded as weakly 
coupled semi-independent degrees of freedom. We will 
discuss this further in Sec. IIXI 

For systems very close to the percolation point, p* ~ 
0.407, one cannot expect to detect differences from the 
behavior exactly at p* . For the smaller cluster sizes at 
p = 0.3 we can observe in Fig. [M] what is likely a cross- 
over behavior from the behavior at p* to the asymptotic 
scaling behavior at p = 0.3. The effective exponent below 
sizes L « 20 is smaller than the value we found at p*, 
but, on the other hand, p = 0.3 is already quite far away 
from p* and it is not surprising that a different behavior 
obtains here. Closer to p* we expect data for small sizes 
to scale with z « 2D/, but to observe clearly this scaling, 
followed by a cross-over to z = D — 2, would require 
larger clusters than we can access currently. 

It should be noted that the results discussed here (and 
those for the bilayer in the previous section) do nei- 
ther prove that the mapping to quantum rotors holds for 



p < p* (and in the bilayer at p*), nor that the dynamic 
exponent exactly equals D = 2 (or Df). For fracton ex- 
citations, one would expect z ^ Df (but close to Df);^ 
It would therefore be useful to determine z for the single 
layer at p < p* and the bilayer at p* to higher precision, 
which, however, is a very demanding task that we leave 
for future studies. 



VII. CLASSICAL DIMER-MONOMER MODEL 

In the mapping of a quantum antiferromagnet onto a 
quantum rotor model, one assumes that there is local 
antiferromagnetic order on some length scale A. A sub- 
system i of the system, of length A, is then replaced by 
a quantum rotor L^, which can reproduce the "Ander- 
son tower" of low-energy states of different total spin S 
(which the subsystem would exhibit in isolation). The 
rotors for all the subsystems are then coupled in a way 
consistent with the expected dominant fluctuations and 
symmetries of the system. For such a mapping to pro- 
duce the correct physics, the subsystems should consist 
of an even number of spins, arranged in such a way that 
their ground state, in isolation, is a singlet. If the sys- 
tem geometry does not allow for such a decomposition, 
the situation will be more complicated. The question is 
then; how can one decompose the system into quantum 
rotors and "left-over" spins in a well defined way, which 
maintains the salient features of the disordered clusters? 

The smallest unit for which a local quantum rotor can 
be considered in a hypothetical mapping is a dimer con- 
sisting of two nearest-neighbor spins, which in isolation 
has a singlet (/ = 0) ground state and a triplet {I = 1) 
excited state. This corresponds to a quite severely trun- 
cated rotor tower, but the local cut-off should not mat- 
ter for the low-energy physics of the coupled system. We 
have already discussed the fact that a disordered cluster 
cannot normally be fully decomposed into such dimers, 
as there would in most cases be some "dangling spins" 
(or, more generally, regions of imbalance in the sublattice 
occupation numbers) left over after the cluster has been 
maximally covered with close-packed dimers. If we con- 
sider larger subsystems, there will be similar problems, 
i.e., not all subsystems will have singlet ground states 
in isolation. We will here proceed to investigate the geo- 
metric decomposition of the system into nearest-neighbor 
dimers and left-over monomers. 

In a standard classical dimer model;^ a dimer cor- 
responds to two connected nearest-neighbor sites, here 
on a square lattice. The statistical mechanics problems 
corresponds to counting all the dimer coverings. In a 
dimer-monomer model, '^^ there are also some unpaired 
sites present, and the counting now includes all possible 
dimer and monomer configurations; normally at a fixed 
density of monomers. In the case at hand here, we in- 
vestigate disordered clusters, and we want to maximally 
cover the clusters with dimers. There will then typically 
be some left over monomers that cannot be paired. For a 
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FIG. 15: (Color online) Updating processes used in Monte 
Carlo sampling of the classical dimer-monomer aggregation 
model, (a) is the annihilation of two monomers, leading to a 
dimer. (b) shows the two elementary monomer-dimer moves. 
In (c), a monomer pair is temporarily created out of a dimer. 
One of the monomers is then moved until it can be annihilated 
with another monomer. A large number of dimers can be 
changed in such "loop updates". 



given cluster, we want to sample dimer-monomer config- 
m'ations with the smallest possible number of monomers. 
We are interested in the spatial distribution of monomers, 
which provides us with a concrete quantitative measure of 
"sublattice imbalance" . We want to identify the regions 
of sublattice imbalance and investigate the size distribu- 
tion of these regions. 

Here we consider clusters constructed on L x L lat- 
tices, with, as before, only the largest cluster found in 
each realization included in the statistics. In Monte Carlo 
sampling of the dimer-monomer configurations on these 
clusters, we start with an arbitrary configuration, e.g., 
one containing only monomers. The updating scheme is 
illustrated in Fig. 1151 In (a), when two monomers are lo- 
cated next to each other, they are annihilated and form 
a dimer. Dimers and monomers can be updated together 
according to the two moves shown in (b). We can also 
break a dimer into two monomers, as in (c). One of the 
monomers is then moved, together with dimers as in (b), 
until it encounters a monomer (which can, but does not 
have to be, the same as the one it was originally paired 
with), together with which it again can be combined to 
form a dimer as in (a). This is an efficient way to up- 
date parts of the cluster where there are no monomers 
(other than the two introduced for the purpose of the up- 
date). This simulation process will eventually converge 
to a state with the minimum number of monomers for a 
given cluster, because the monomer annihilation process 
(a) is always carried out when possible. Whenever two 
monomers are created, they will eventually be annihi- 
lated. Our model is therefore also an aggregation model 
for dimers. 

For a given cluster generated on an L x L lattice, af- 
ter a long equilibration to make sure that the minimum 
monomer number has been reached, we collect statistics. 
One quantity of interest is the average monomer density 
for each site. Most of the sites never have any monomers. 
We define a "moment" as a region consisting of Sm sites 




FIG. 16: (Color online) Probability distribution of the mo- 
ment size of the classical dimer-monomer model at the per- 
colation point p* (top panel) and at p — 0.3 (bottom panel) 
graphed on a log-log scale. Hundreds of realizations of the 
largest cluster on L x L lattices were used. For large L the 
distributions collapse onto a single curve, reflecting a finite 
typical moment size. The average moment size is (Sm) = 16 
at p* and 47 at p = 0.3. The asymptotic form of the proba- 
bility distribution for large clusters is oc e"^""''^, as shown in 
the insets using a semi-logarithmic scale. 



to which one or several monomers are confined. By defi- 
nition, a monomer inside such a moment cannot move to 
a different moment through the Monte Carlo processes. 
In addition, the moments consist only of sites on the same 
sublattice, because an individual monomer only moves on 
a given sublattice, as in Fig. llSr b). Two moments on dif- 
ferent sublattices cannot have any sites that are nearest 
neighbors, because then two monomers in these different 
regions could become adjacent and annihilate each other. 

Keeping track of the moment regions and their sizes 
involves straight-forward book-keeping, and we just pro- 
ceed to discuss results. Fig. [T6l shows the size distribution 
of the moments both at the percolation point and away 
from it, at p = 0.3. For small moment sizes Sm, the dis- 
tribution is close to a power-law, especially at p — 0.3, 
but there is a cross-over to a clearly exponential decay 
for large Sm ■ The distribution can be fitted well with the 
form e-'^'"/'^, with ct « 42 and « 300 at p* and p = 0.3, 
respectively. The average moment size (Sm) computed 
as a sum over all the sizes is smaller; {Sm) ~ 16 at p* 
and « 47 at p = 0.3. In the figure, the largest cluster 
sizes are much larger than (Sm) and cr, and the curves 
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FIG. 17: (Color online) (a) Examples of the convergence of 
the average triplet IPR for three individual clusters as a func- 
tion of the projection power, here normalized by the cluster 
size as P/n. The clusters were generated on 26 x 26 lattices. 
The dashed curves are of the form a + be~''^^", with a, b, c 
adjusted to fit the last few (large P/n) points. 



for the largest L overlap almost completely. 

These calculations prove that the notion of local sub- 
lattice imbalance is well defined and quantifiable. Finite 
moment regions exist both at and away from the perco- 
lation point. In the next section, we will present result 
from valence bond quantum Monte Carlo simulations in 
the triplet sector. We will there also look at the spa- 
tial distribution of the monomers in the classical dimer- 
monomer model, and compare it with the distribution of 
the triplet in the lowest excitation of the actual quantum 
spin model. 



VIII. SPATIAL DISTRIBUTION OF 
SINGLET-TRIPLET EXCITATIONS 



As discussed in sec. Ill C[ the valence-bond projector 
QMC method offers us the possibility to examine the 
lowest triplet state in a unique way. In a disordered sys- 
tem, the spatial distribution of the triplet bond gives a 
very direct measure of the extent to which different parts 
of the system are affected when exciting a cluster with 
S = ground state to its lowest S* = 1 state. An example 
of the triplet density for a very small cluster was already 
presented in Fig.[TJ It should be noted that the statistics 
of the singlet bonds is also affected by the presence of a 
triplet bond, and, thus, just examining the properties of 
the triplet bond does not give a complete picture of the 
excitation. However, if a large region of the system has 
no (or very low) average triplet density, then the singlets 
of that region should also not be much affected. The spa- 
tial distribution of the triplet should therefore provide a 
valid measure of the tendency (if any) of the triplet ex- 
citation to localize. 

We study the site dependent triplet density pi = 
{nt{i)), where the triplet occupation number nt{i) is de- 
fined such that if a triplet bond connects sites i and j. 




(c) 



£ 



FIG. 18: (Color online) Properties of a typical cluster with 
288 sites; the classical monomer density (a), the triplet den- 
sity (b), and the local gap estimate Ai (c). All quantities are 
shown on a color scale ranging from the smallest (0 in the 
case of pt and pm) to the largest value. The absolute values 
are irrelevant for the purpose of the discussion here. In (a), 
the numbers inside the squares label the different classical 
moment regions. The dots indicate empty sites. 



then nt{i) = 1 and nt(j) = 1, while nt{k) = for all other 
sites k. In addition to visually examining the triplet den- 
sity for representative individual clusters, it is also useful 
to have a quantitative measure of localization. For this 
purpose, we use the inverse participation ratio (IPR) cor- 
responding to the triplet density; 



Rt = 



En 



En 



^=lPi 



(42) 
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FIG. 19: Log-log plot showing the finite-size scaling of the 
IPR of the triplets {Rt) and the classical monomers (Re)- 
The lines correspond to scaling L? , with 7 = 1.39(3) (based 
on a line fit) for Rt and 7 = Dj for R^. 



This quantity characterizes the number of sites involved 
in a triplet excitation, and can be averaged over cluster 
realizations. Two extreme cases can help to clarify the 
meaning of i?^; if the triplet is completely localized on 
only two sites, then Rt = 2, while if it is equally spread 
out over all the sites of an n-site cluster, then Rt — n. 
We will study the dependence of (Rt) on the cluster size. 

In a projector method based on a power , one con- 
verges to the lowest state in a given symmetry sector 
when the power P of the hamiltonian is sufficiently high. 
For an n-site cluster, one would expect that the P re- 
quired for convergence scales as n or worse. This can 
be seen if we compare with an alternative projection 
method — the imaginary-time evolution e~''^|4') of the 
trial state. Here /? is analogous to an inverse temper- 
ature; starting from a state at some "temperature", we 
"cool it" by increasing /3. A better trial state corresponds 
to a lower initial temperature. For large /?, the dominant 
power P in a Taylor expansion of e~^^ is P — P\Eq\, 
where Eq is the ground state energy, which is propor- 
tional to the cluster size n. Thus, if we project with just 
a fixed power P oi H, we would get essentially the same 
result if P « P (X /3n. The energy scale of the excitations 
decrease with increasing n (very quickly so in the problem 
under consideration here, because the dynamic exponent 
is large), and we should therefore expect to need larger 
f3 for larger n. Thus, in the fixed-power scheme, the P 
required for convergence should increase as some power 
(larger than one) of n. 

We show examples of the convergence of the IPR for 
three different clusters in Fig. [T71 The large fluctua- 
tions in the finite-size gap, discussed in Sec. lIIIl naturally 
also imply large variations in the convergence rate of the 
triplet IPR (which is governed by the gap between the 
first and second triplet, which also exhibits large fluctu- 
ations). As explained in Sec. IlICi we are restricted to 
P for which the triplet survival probability in the pro- 
jection is reasonably large. In order to ensure that the 
results truly reflect the lowest excitation for each clus- 
ter, we carry out extrapolations to infinite P using a 
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FIG. 20: (Color online) Probability distribution of the length 
{dx,dy) of the singlet (top panel) and triplet (bottom panel) 
bonds for the cluster in Fig. 1181 The singlet distribution is 
strongly peaked at short bonds, and we have therefore cut off 
the corresponding peak in the low left corner of the histogram. 
The the remaining weight is 9% of the total. 



simple exponential form, as explained in the caption of 
Fig. [T71 The fluctuations of the disorder averaged (Rt) 
are completely dominated by the cluster-to-cluster vari- 
ations, and we believe that any remaining errors related 
to the convergence are smaller than the final error bars 
(based on a few hundred clusters of each size). 

We first examine the spatial distribution of the triplets. 
The triplet density for each site of a typical cluster is 
shown using a color scale in Fig. llSf b). Here we compare 
the triplet density with two other calculations — the clas- 
sical monomer density in (a) and the local gap A.^ in (c). 
It is apparent that the triplet is concentrated to a rela- 
tively small fraction of all the sites of the cluster. At the 
same time, the affected sites form groups that are spread 
out over the cluster. This is exactly in agreement with 
our hypothesis of low-energy excitations involving a num- 
ber of localized moments. It is also clear from Fig.fTSlthat 
the classical monomer density is high wherever the triplet 
probability is significant. This proves that our measure 
of sublattice imbalance in terms of classical monomers 
indeed corresponds very closely to the actual locations 
affected by excitations. Note also that some sites with 
high monomer density do not have a high triplet den- 
sity. This is also expected, because the lowest triplet 
excitation should not necessarily involve all of the classi- 
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cal monomer regions. Higher triplet states may involve 
other subsets of moments. Finally, there is a very good 
correspondence between regions of low local gaps and 
hight triplet density. 

Next, we discuss the IPR of the triplet. It is interest- 
ing to compare this with the total number of spins in the 
classical dimer-monomer model. We therefore also define 
a classical IPR, as in Eq. (|42|) but with the triplet density 
Pi replaced by the classical monomer density. Both these 
IPRs, averaged over several hundred clusters, are shown 
versus the cluster length L on a log-log scale in Fig. [191 
They both scale according to power laws. The classical 
IPR is consistent with the form L^^ ~ (n). In combi- 
nation with the fact that the individual moment regions 
are finite, as we showed in the previous section, this is in 
agreement with our extremal- value analysis in Sec. IIIIB| 
which relied on the number of effective moments being 
proportional to n. However, the triplet IPR scales with 
a smaller power; ~ with 7 = 1.39(3). Thus, 

not all the effective moments are involved in the lowest 
excitation, but since the size of the excitation still grows 
with L these are not localized excitations. 

Another important aspect of the valence-bond calcu- 
lation is that the bond lengths also contain information 
directly pertaining to the nature of the excitations. In 
Fig. [20] we show the distributions of both the singlet and 
triplet bonds lengths for the cluster shown in Fig. [TH| 
While the triplet bond is typically long, the singlet dis- 
tribution is strongly peaked for the shortest bonds. We 
have therefore cut off more than 90% of the weight in 
the singlet histogram in order to be able to show the 
more interesting distribution of long bonds. Every peak 
in the triplet distribution can be perfectly matched to 
the distance between two regions (on different sublat- 
tices) with a high concentration of triplets/monomers in 
Fig. [T8| This again supports the notion of excitations 
of weakly interacting effective moments. In the singlet 
distribution, there are also features corresponding to the 
same lengths as in the triplet case. This is also what one 
would expect if the triplet state essentially corresponds 
to promoting a long singlet to a triplet in a superposi- 
tion, as discussed in Sec. lIVCl The average length of the 
triplet bond also scales with the cluster size according 
to a power-law, as shown in Fig. 1211 This power law, in 
combination with that for the triplet IPR (Fig. [TO]) and 
classical percolation exponents, should be related to the 
dynamic exponent z « 3.6. Exactly how is not presently 
clear, however. 



IX. SUMMARY AND DISCUSSION 

To summarize, we have discussed several calculations 
aimed at elucidating the quantum dynamics of the S ~ 
1/2 antiferromagnetic Heisenberg model on randomly di- 
luted clusters. Quantum Monte Carlo simulations in 
combination with sum rules show that the low-energy ex- 
citations at the percolation point are described by an un- 
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FIG. 21: (Color online) Size dependence of the average length 
of the triplet bond on L x L clusters. The line corresponds to 
a power-law divergence Lt ^ L" , with a — 0.81 ± 0.03. 



expectedly large dynamic exponent; z « 3.6±0.1 « 2D/, 
where Df = 91/48 is the fractal dimension of 2D perco- 
lation. Using extremal-value statistics, we were able to 
relate z and two exponents characterizing the probabil- 
ity distribution of local gaps (an exponent a governing 
the size dependence and uj describing the distribution of 
small gaps) , according to Eq. ([55]) . This kind of scaling 
indicates that the excitations involve effective localized 
finite magnetic moments, which interact through the re- 
maining, magnetically inert parts of the cluster. This is 
also confirmed directly by imaging the spatial distribu- 
tion of the triplet excitations in the valence bond basis, 
where the triplet state can be described in terms of a pair 
of spins forming a triplet in a "singlet soup" of valance 
bonds. The triplet bond fluctuates between several iso- 
lated regions, and its average length scales as a power of 
the cluster size. The average number of spins affected 
by the excitation also grows as a power of the system 
size, with (3 ~ 0.74. The excitations are thus localized at 
multiple moment regions, which are spread out over the 
cluster. All these results lead to a picture of an effective 
low-energy system consisting of a network of globally en- 
tangled local moments, where the moments correspond 
to regions of sublattice imbalance. 

We have introduced a quantitative measure of sublat- 
tice imbalance, in terms of a classical dimer-monomer ag- 
gregation model. Monte Carlo simulations of this model 
show that the monomers form isolated finite regions, and 
the number of such regions scales linearly in the cluster 
size n. Sites with a high triplet concentration coincide 
very well with high monomer density, confirming directly 
that sublattice imbalance in the Heisenberg model is as- 
sociated with the formation of weakly interacting effec- 
tive moments. 

We have also shown that when two identical clusters 
are coupled in a bilayer, with a small inter-layer cou- 
pling J± (smaller than the value at which the long-range 
order vanishes^), the low-energy excitations change dra- 
matically. A finite-size scaling analysis for clusters with 
J_\_/J = 0.01 show a much smaller dynamic exponent, 
z K Df, than for the single-layer clusters. There is 
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no sublattice imbalance in this "dimer diluted" bilayer 
model, and the results therefore provide additional evi- 
dence for the important role played by effective moments 
at imbalanced regions in the single-layer clusters. 

The result z « 3.6 was obtained by studying clusters 
in which the sublattices are balanced globally, i.e., the 
ground state spin is S* = 0. We have also pointed out 
that for clusters with global imbalance in the sublattice 
occupation (leading to 5 > 0), the dynamic structure 
factor has spectral weight predominantly arising from the 
S ^ S ~\- \ channel, while there is very little weight in 
the S ^ S and S ^ S — 1 channels (apart from the 
elastic 5* — > 5* weight). Experiments directly probing 
the inelastic spectral weight, e.g., neutron scattering and 
nuclear magnetic resonance, should be dominated by the 
S ^ S + 1 channel, and S should be typically large (oc 
y/n) for random clusters. For this situation our sum rule 
method gives a smaller effective dynamic exponent, z « 
1.5Df, than for the = — > 1 excitations of globally 
balanced clusters. The lowest-energy excitations, which 
are in the S ^ S —1 channel, are not accessible with the 
sum rule approach. We have argued, based on an analysis 
of approximate (variational) valence bond states, that 
their energy should scale with the same z w 3.6 ~ 2Df 
as in the case of = ^ 1 excitations. 

There have been attempts previously to determine the 
dynamic exponent of the diluted Heisenberg model based 
on quantum Monte Carlo calculations. Yu et al. studied 
the temperature dependence of the correlation length and 
concluded that it scales in a way consistent with z = Df 
at the percolation point However, the scaling assump- 
tion was one corresponding to quantum-criticality, which 
may not apply because the percolating cluster at T = 
does not have quantum critical fluctuations in the sense 
of power-law decaying correlation functions. A similar 
treatment of the clean 2D Heisenberg model would fail 
to give z = D (which is not a quantum-critical exponent 
but one characteristic of the quantum-rotor excitations 
of the Neel state) , because the correlation length diverges 
exponentially as T ^ Oi^ 

We have here focused exclusively on the dynamics of 
the percolating cluster. In order to relate the results 
quantitatively to specific experiments, one should include 
the contributions from all clusters. The cluster distribu- 
tion is given by classical percolation theory^i^ which can 
be combined with the finite-size scaling properties that 
we have found here for the distributions of the local and 
smallest gaps. We have also not discussed the conse- 
quences of our r = results for the T > behavior. 
This is of course also an important experimental issue, 
and will be interesting to consider in future studies. 

Moving away from the percolation point p*, our cal- 
culations show that the dynamic exponent of the sin- 
gle layer is z Ri 2 = D. However, the classical dimer- 
monomer model has finite localized monomer regions also 
away from the percolation point. This suggests that the 
change in the spin dynamics upon moving away from p* is 
related to the effective interactions between the moment 



regions, not the disappearance of the moments. Such a 
qualitative change in the interaction aspects of the mo- 
ments is not completely unexpected, since the spin stiff- 
ness of the percolating cluster at p* is strictly zero in the 
thermodynamic limit^ i^^i'^^ (although the cluster is or- 
dered), whereas it becomes finite (according to a power- 
law) away from the percolation point. The more robust 
cluster order for p < p* should qualitatively change the 
effective interactions between distant monomer regions, 
likely locking all of them to the global Neel vector (which 
is the case for a single moment in a 2D systemiiii^) . The 
effectively independent nature of the magnetic moments 
exactly at the percolation point (and the very weak in- 
teractions between them, are thus intimately related to 
the fractal structure and related vanishing spin stiffness 
of the network connecting the moment regions. We pre- 
sented some results showing the cross-over from scaling 
controlled by the percolation point to 2D behavior. 

In spite of the close agreement with the dynamic expo- 
nent expected based on the quantum rotor mechanism^ 
for the single layer away from p* and the bilayer at p* , it 
is still not certain that the lowest-energy excitations in 
these systems are quantum rotor states. Over the years 
there have been considerable efforts to understand the 
dynamics of various randomly diluted systems close to 
and at the percolation point. The "fracton" has been 
introduced as a generic excitation which develops out of 
plane waves (e.g., spin waves for an antiferromagnet) for 
a translationally invariant system upon dilutioupi^ Nu- 
merical calculations based on spinwave theory show that 
the dynamic exponents for fractons in the 2D percolating 
antiferromagnet is very close to Dj^^ This calculation 
does not properly account for the vanishing spin stiff- 
ness of the percolating cluster at p* and the existence 
of localized moments, but it may still be valid close to 
the percolation point (where a finite stiffness develops). 
It is possible that the dynamic exponents z D — 2 
and z K, Df that we have obtained here for the single 
layer with p < p* and the bilayer at p*, respectively, 
are due to fractons, not quantum rotors. However, ex- 
actly at the percolation point, the physics of the globally 
entangled local moments that we have discussed here is 
clearly different from fractons (which can exist in sys- 
tems that do not have any objects corresponding to lo- 
calized moments) and the value of the dynamic exponent 
is twice that expected based on fractonSfiS It is thus pos- 
sible that several types of excitations co-exist at and close 
to the percolation point; quantum rotors, fractons, and 
the globally entangled moment excitations. 

One may still be able to describe the low-energy 
physics of the system away from the percolation point 
as a network of weakly interacting moments, but now in 
the presence of a staggered field mimicking the coupling 
to a common Neel vector (i.e., the sign of the field de- 
pends on which sublattice a moment is associated with on 
the original lattice) . The strength of the effective stag- 
gered field (which for a finite cluster should be allowed 
to have a fiuctuating direction as well^*) should increase 
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upon moving away from the percolation point (being at 
p*, due to the vanishing spin stiffness — the energy scale 
of twisting the Neel order globally) . Most likely, even an 
infinitesimal field will asymptotically (for large clusters) 
change the dynamic exponent. 

The role of effectively isolated spins in the formation 
of long-range order= on the percolating cluster has been 
pointed out by Bray-Ali et al.'^^''^^ Although some spins 
can be very weakly coupled (effectively) to the rest of the 
cluster, correlations between them can be stronger than 
within the backbone of the cluster. Arbitrarily weakly 
coupled moments formed by groups of spins can also cor- 
relate over long distances, and hence even a "floppy" frac- 
tal cluster (one with vanishing spin-stiffness) can order 
at T = 0. This picture seemingly contains some of the 
ingredients of our entangled moments picture. However, 
the same ordering mechanism was argued to apply both 
to antiferromagnets and systems of coupled quantum ro- 
tors, whereas we have shown here that the bilayer (which 
should correspond more closely the coupled rotor system, 
since there is no sublattice imbalance) and the single layer 
behave dramatically different. Since the excitations of 
the single layer cluster are also much lower in energy than 
the quantum rotor states considered in previous discus- 
sions of the dynamic exponent,— a theoretical treatment 
within a quantum rotor picture is clearly not adequate. 
In Refs. i38i and i39, the excitations of the diluted system 
were analyzed using spin-wave theory, but this method 
also does not capture the significance of almost isolated 
moments and their long-range global entanglement, and 
no unusually low energy scale was discussed. 

In field-theory language, the "dangling" spins, or re- 
gions of sublattice imbalance, that we have discussed here 
correspond to uncompensated Berry phasesiii Although 
it is quite clear that these should exist in diluted quan- 
tum antiferromagnets, how to properly take them into 
account in analytical calculations for these systems is not 
well understood. To our knowledge, the resulting globally 
entangled moments excitations that we have argued for 
here have not been discussed previously in the literature. 
The effective low-energy system is similar to the random 
antiferromagnet considered by Bhatt and Lee^^ and also 
by Sachdev and Ye4i However, there is an important dif- 
ference in that the nearest-neighbor interactions in our 
system are not frustrated. An effective low-energy hamil- 
tonian should then also not be frustrated. 

The Bhatt-Lee calculation^ was focused on the ther- 
modynamic properties and did not address the dynamic 
exponent. The method applied was a generalization of 
the strong disorder renormalization (singlet decimation) 
scheme by Ma, Dasgupta, and Hu,^^ which has been 
applied to numerous random antiferromagnetic Heisen- 
berg systems j^^i^i^ It would be interesting to apply this 
method also to the diluted clusters. However, there is a 
technical problem in doing this directly, since the decima- 
tion scheme is based on random couplings (successively 
eliminating the strongest coupled spin pair and including 
their remaining effects as modified couplings calculated 



perturbatively) , whereas in the diluted system all cou- 
plings are the same. It may be possible to carry out a 
decimation procedure by eliminating strongly correlated 
spins, instead of strongly coupled ones. The correlations 
could be computed perturbatively based on regions of a 
small number of spins, or using quantum Monte Carlo 
simulations. This way, one could study the renormal- 
ization flows of the correlations and how they relate to 
the sublattice imbalance that we have quantified here in 
terms of the classical dimer-monomer systems. The final 
stages of the decimation procedure should lead to bonds 
(entanglement) between the sites on which our projector 
QMC calculations give a high triplet probability. How- 
ever, we have shown that there are large fluctuations in 
the long bonds (singlet as well as triplet) and it is there- 
fore clear that the scheme cannot asymptotically give the 
correct ground state and low-energy excitation in terms 
of a single bond configurations. In one dimension, the fi- 
nal "random singlet state" is known to be asymptotically 
exact)^ in the sense that a single bond configuration is 
a good representation of a superposition including fluc- 
tuations around this reference state4^ With large fluc- 
tuations of long valence bonds among many moments in 
the percolating clusters, it seems unlikely that a single 
reference configuration would be a good approximation 
in this case. It would still be interesting to investigate 
the flow of the renormalizcd coupling distribution. 

To go further in developing an understanding of the ex- 
citations of the weakly interacting effective moments, in- 
stead of working with the full percolating clusters it may 
be better to explicitly construct the effective low-energy 
hamiltonian we have discussed here. While the geometri- 
cal locations of the moments could be obtained using the 
classical dimer-monomer model, the effective couplings 
are more challenging. One approach would be to just 
study a bipartite network of spins with some suitable 
form of the interactions (which should be non-frustrated, 
with antiferromagnetic couplings between sublattices and 
ferromagnetic intra-sublattice couplings). In principle 
the spins should have mixed S. In one dimension such 
a system is known to have different properties than 
the random 5=1/2 chain with only antiferromagnetic 
couplings.^'' The effective system could be studied with 
the methods used here, as well as with the strong-disorder 
decimation scheme. Comparing results for the moment 
network with the dynamic exponents we have extracted 
here for the full cluster system (and investigating the ro- 
bustness of the exponents to variations in the couplings) 
could shed further light on this challenging problem. 
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